Influence of circadian clocks on adaptive immunity and vaccination responses

The adaptive immune response is under circadian control, yet, why adaptive immune reactions continue to exhibit circadian changes over long periods of time is unknown. Using a combination of experimental and mathematical modeling approaches, we show here that dendritic cells migrate from the skin to the draining lymph node in a time-of-day-dependent manner, which provides an enhanced likelihood for functional interactions with T cells. Rhythmic expression of TNF in the draining lymph node enhances BMAL1-controlled ICAM-1 expression in high endothelial venules, resulting in lymphocyte infiltration and lymph node expansion. Lymph node cellularity continues to be different for weeks after the initial time-of-day-dependent challenge, which governs the immune response to vaccinations directed against Hepatitis A virus as well as SARS-CoV-2. In this work, we present a mechanistic understanding of the time-of-day dependent development and maintenance of an adaptive immune response, providing a strategy for using time-of-day to optimize vaccination regimes.


Rhythmic homing dynamics of T cells (ID1)
The mathematical model used to describe the homing and egress dynamics of T cells for skin-specific lymph nodes is based on the model described in Druzd el al. 7 . In this previous study, it was found that homing and egress of lymphocytes are both assumed to follow sinusoidal oscillating dynamics. The lymphocyte count within an individual lymph node, X(t), therefore changes over time according to Hereby, A and d define the homing and egress rates of cells, respectively, and ω the angular frequency of the oscillating dynamics given a 24h-cycle. The individual phases of the time-dependent gain and loss rates of lymphocytes are shifted by the parameters ϕ h and ϕ e , respectively.
The model in Eq.
(1) was fitted to the experimental data on the lymphocyte counts of the axillary, inguinal and superficial cervical lymph node 7 using the deSolve and optim-package in R. Parameter estimates for the individual lymph nodes are given in Table S1, with 95%-confidence intervals of estimates obtained by profile likelihood analysis 33 . Individual data and best model predictions are shown.

Homing dynamics of Dendritic cells (ID2)
The homing dynamics of dendritic cells (DC) was determined using the data of the crawl-In assay given in Holtkamp et al. 18 . In brief, labelled bone marrow-derived DC were injected into ears of mice that were harvested at different ZT times. The infiltration of DC into the vasculature was then followed by 2-photon microscopy. Data indicate a ZT-dependent infiltration of cells into the vasculature. Distinguishing between the extra-and intravascular concentration of DC, i.e., DC E and DC I , respectively, the dynamics after injection is then described by the following system of ordinary differential equations: Hereby, the parameter α defines the amplitude of the rhythmic component of the homing rate with angular frequency ω given a 24h-cycle and phase-shift parameter ϕ h,DC , while μ specifies a constant homing rate independent of the time of injection and ear harvest.
The model given by Eqs. (2) and (3) was fitted to the observed ratio of intra-and extravascular cells using the deSolve and optim-package in R. The initial concentration of DC injected was set to 1 without loss of generality, DC E (0) = 1. Parameter estimates are given in Table S1, with 95%-confidence intervals of estimates obtained by profile likelihood analysis 33 . Individual data and best model predictions are shown.

T cell proliferation dynamics (ID3)
To determine the observed rhythmic within the proliferation of T cells (Figure 3c, main manuscript), we describe the time-dependent oscillating dynamics of the division index, DI, relative to ZT7 by Hereby, DI 7 defines the baseline division index at ZT7, and Λ, ω, and ϕ P the amplitude, angular frequency and phase-shift of the rhythmic proliferation component, respectively. Parameters were estimated using a maximum likelihood approach as before in R, with 95%-confidence intervals of estimates obtained by profile likelihood analysis 33 . Data and best model predictions are shown with parameter estimates given in Table S1.

Mathematical model describing interaction of multiple rhythms (ID4)
With rhythmicity detected in the individual components of lymphocyte and dendritic cell (DC) homing to individual lymph nodes (LNs), as well as in cell proliferation and time-dependent cell velocity dynamics within the LN, we developed a mathematical model to determine how the interaction of the different rhythmic components shapes immune responses and leads to maintained rhythmicity. To this end, we combined the individual dynamics shown in Eqs.
(1) -(4) before. We additionally assume that a concentration of activated T cells, T A , is formed based on the interaction of lymphocytes and DC in the LN at an activation rate σ. This activation is shaped by a factor ν ZT , which accounts for the different velocities of T cells dependent on ZT with ν ZT =1 for ZT=7. Activated T cells start to proliferate, T P , at rate ρ, which is additionally shaped by the rhythmic division index, DI(t) as defined within Eq.
(4). To prevent unlimited LN expansion, a carrying capacity N is assumed that defines the maximal number of cells in a lymph node. Therefore, infiltration of unactivated T cells, X(t), as well as proliferation of activated T cells is assumed to saturate with increasing LN expansion. In addition, DC are lost from the LN at a rate γ. A sketch of the model is provided in Figure 4b of the main manuscript. The whole model is then defined by: d (1 + sin(!(t + e ))) In comparison to the description of lymphocyte homing and egress in Eq.(1), we found that two modifications for the model were necessary to describe the observed LN expansion, as well as maintaining rhythmicity throughout time. These modifications account for altered lymphocyte recruitment and egress dynamics due to T-DC interactions within the lymph node. The function f 1 (ω,ϕ h ,t) accounts for the fact that T-DC interactions facilitate lymphocyte recruitment (e.g. due to inflammation) by ensuring lymphocyte increase to be not smaller than the influx generated at the time of injection of DC, i.e., t ZT . Furthermore, a delayed egress of activated lymphocytes within in the lymph node is assumed, with the egress rate d reduced in a density-dependent manner. Hereby, the parameter λ denotes a scaling factor.
The model was parameterized using the estimates of the individual components obtained from ID1-ID3 (see Table S1). In addition, the velocity factor ν ZT was calculated based on the ratio of the mean velocities measured (Figure 1f). Furthermore, we assumed an activation rate of σ = 0.011 h -1 and a cell proliferation rate of ρ = 0.03 h -1 , corresponding to an average division period of~33 h, to match reasonable cell numbers within LN. The average dwell time of DC within the LN contributing to T cell activation, 1/γ, was set to 50 h, and the scaling factor regulating the efflux of activated T cells was set to λ = 10 -3 based on the analyses of previous data sets. All individual parameters used for simulation are given in Table S1.
The different skin draining lymph nodes, i.e. axillary, inguinal and superficial cervical LN, were simulated individually and the average dynamics across all lymph nodes was calculated to compare LN expansion using different ZT time points for DC injection. DC were assumed to migrate to the LN for a time period of 1h after injection. We tested different model assumptions with the mentioned feedback dynamics of T-DC interactions on lymphocyte recruitment and egress being able to explain the observed maintenance of rhythmicity within the immune response at later time points.